Загрузите данные в датафрейм.

gmp <- read.table("https://raw.githubusercontent.com/SergeyMirvoda/MD-DA-2017/master/data/gmp.dat")
gmp$pop <- gmp$gmp/gmp$pcgmp
estimate.scaling.exponent <- function(a, y0=6611, response=gmp$pcgmp,
                                        predictor = gmp$pop, maximum.iterations=100, deriv.step = 1/100,
                                        step.scale = 1e-12, stopping.deriv = 1/100) {
  mse <- function(a) { mean((response - y0*predictor^a)^2) }
  for (iteration in 1:maximum.iterations) {
    deriv <- (mse(a+deriv.step) - mse(a))/deriv.step
    a <- a - step.scale*deriv
    if (abs(deriv) <= stopping.deriv) { break() }
  }
  fit <- list(a=a,iterations=iteration,
              converged=(iteration < maximum.iterations))
  return(fit)
}
a<-0.15
expr <- estimate.scaling.exponent(a)
estimate.scaling.exponent(a)
## $a
## [1] 0.1211533
## 
## $iterations
## [1] 58
## 
## $converged
## [1] TRUE

С помошью полученного коэффициента постройте кривую (функция curve) зависимости

plot(pcgmp~pop, data=gmp, xlab="Население", log="xy",ylab="Доход на душу населения ($/человеко-год)", main="Метрополии США, 2006")
curve(6611*x**expr[[1]],add=T,col='red')

Удалите точку из набора исходных данных случайным образом, как изменилось статистическая оценка коэффициента a?

deleting_row <- function(a,del.row = 0){
  if(del.row==1){
    row.for.del <- sample(nrow(gmp),1)
    K <-  which(!(1:nrow(gmp)) %in% row.for.del)
    gmp <<-gmp[K,]
  }
  gmp$predicted <- 6611*(gmp$pop)^estimate.scaling.exponent(a)[[1]]
  message(paste("a:",estimate.scaling.exponent(a)[[1]]," iterations:", estimate.scaling.exponent(a)[[2]],"\n"))
  lines(gmp$pop, gmp$predicted, col = "red")
  #return(estimate.scaling.exponent(a)[[1]])
}

Запустите оценку несколько раз с разных стартовых точек. Как изменилось значение a?

o1<-c()
i=1
for (u in seq(0.20,0.4,by = 0.01)){
  o1[i]<-estimate.scaling.exponent(u)[[1]]
  i<-i+1
  plot(pcgmp~pop, data=gmp, xlab="Население", log="xy",ylab="Доход на душу населения ($/человеко-год)", main="Метрополии США, 2006")
  deleting_row(u)
}
## a: 0.121153292739596  iterations: 70

## a: 0.121153292739592  iterations: 78

## a: 0.121153292739593  iterations: 100

## a: -0.0158227947406229  iterations: 100

## a: -0.282180595453594  iterations: 100

## a: -0.475270609311027  iterations: 100

## a: -0.719596634464836  iterations: 100

## a: -1.04635641522945  iterations: 100

## a: -1.48417713043446  iterations: 100

## a: -2.06949635226438  iterations: 100

## a: -2.85063412085794  iterations: 2

## a: -3.89182201309161  iterations: 2

## a: -5.27849638056654  iterations: 2

## a: -7.12437014241724  iterations: 2

## a: -9.58089019392308  iterations: 2

## a: -12.8498946401176  iterations: 2

## a: -17.2005650956782  iterations: 2

## a: -22.9921483818708  iterations: 2

## a: -30.7044336419058  iterations: 2

## a: -40.9786620271726  iterations: 2

## a: -54.6724802138931  iterations: 2

plot(o1,main="Оценка а с разных стартовых точек")

for(i in 1:100){
  #print(paste("Количество точек", nrow(gmp)))
  plot(pcgmp~pop, data=gmp, xlab="Население", log="xy",ylab="Доход на душу населения ($/человеко-год)", main="Метрополии США, 2006")
  deleting_row(0.12,1)
}
## a: 0.121114444808736  iterations: 54

## a: 0.121109598973552  iterations: 54

## a: 0.121051791424099  iterations: 54

## a: 0.121063122493425  iterations: 54

## a: 0.121073882368287  iterations: 54

## a: 0.120962639733348  iterations: 54

## a: 0.120997617266159  iterations: 54

## a: 0.120982488200649  iterations: 53

## a: 0.120958073079915  iterations: 54

## a: 0.120969302525322  iterations: 54

## a: 0.120906250499022  iterations: 54

## a: 0.120951431111156  iterations: 54

## a: 0.120853195998464  iterations: 54

## a: 0.120802341916282  iterations: 53

## a: 0.120733852796318  iterations: 53

## a: 0.1204034114657  iterations: 53

## a: 0.120399875214202  iterations: 53

## a: 0.120556015906425  iterations: 53

## a: 0.120564855239912  iterations: 53

## a: 0.120491635720752  iterations: 53

## a: 0.120509338571741  iterations: 53

## a: 0.120504228779825  iterations: 53

## a: 0.120504165847206  iterations: 53

## a: 0.120549051575398  iterations: 53

## a: 0.120493073216661  iterations: 53

## a: 0.120614826927809  iterations: 53

## a: 0.120637954086077  iterations: 53

## a: 0.120611056710325  iterations: 53

## a: 0.120598440844703  iterations: 53

## a: 0.120537423988755  iterations: 53

## a: 0.12046843221171  iterations: 53

## a: 0.120509538771279  iterations: 53

## a: 0.120529594053472  iterations: 53

## a: 0.120508036439313  iterations: 52

## a: 0.120506086702969  iterations: 52

## a: 0.120499567611679  iterations: 52

## a: 0.120493177058957  iterations: 52

## a: 0.120479958735048  iterations: 52

## a: 0.120369508449169  iterations: 52

## a: 0.120386125745304  iterations: 52

## a: 0.120367624500289  iterations: 52

## a: 0.120349299709309  iterations: 52

## a: 0.120387027871802  iterations: 52

## a: 0.120372057746726  iterations: 52

## a: 0.12039504646698  iterations: 52

## a: 0.120407686223157  iterations: 52

## a: 0.120411112022311  iterations: 52

## a: 0.120341609610505  iterations: 51

## a: 0.120281050204919  iterations: 51

## a: 0.12032266622458  iterations: 51

## a: 0.120325902335481  iterations: 51

## a: 0.120379271312484  iterations: 52

## a: 0.120462893258827  iterations: 52

## a: 0.120468834620813  iterations: 52

## a: 0.120516915913398  iterations: 52

## a: 0.120460627673639  iterations: 52

## a: 0.120537524987198  iterations: 52

## a: 0.120516912270523  iterations: 52

## a: 0.12057491253572  iterations: 52

## a: 0.120612389670403  iterations: 52

## a: 0.120586085165935  iterations: 53

## a: 0.120506010814773  iterations: 52

## a: 0.120603503069724  iterations: 53

## a: 0.12057187863211  iterations: 53

## a: 0.120589879005083  iterations: 53

## a: 0.120620591644072  iterations: 53

## a: 0.120651081755946  iterations: 53

## a: 0.120624995371933  iterations: 53

## a: 0.120519862402356  iterations: 53

## a: 0.120558014926904  iterations: 53

## a: 0.120641067954741  iterations: 53

## a: 0.120903908632093  iterations: 54

## a: 0.120839183151217  iterations: 54

## a: 0.120884075745674  iterations: 54

## a: 0.120893377387248  iterations: 54

## a: 0.120881914533995  iterations: 53

## a: 0.120862005846508  iterations: 53

## a: 0.120848245022099  iterations: 53

## a: 0.120918111218219  iterations: 53

## a: 0.120845560403798  iterations: 53

## a: 0.120798376832012  iterations: 53

## a: 0.120638575858814  iterations: 53

## a: 0.120378650874897  iterations: 52

## a: 0.120442199792226  iterations: 53

## a: 0.120401481018942  iterations: 53

## a: 0.120388054749061  iterations: 53

## a: 0.120398773523536  iterations: 53

## a: 0.120349251907562  iterations: 52

## a: 0.120367062479791  iterations: 53

## a: 0.12039077103302  iterations: 53

## a: 0.12042560000916  iterations: 53

## a: 0.120317479135247  iterations: 52

## a: 0.120254098775309  iterations: 52

## a: 0.120209377634074  iterations: 51

## a: 0.120272847012106  iterations: 52

## a: 0.120234784940414  iterations: 51

## a: 0.120250921498287  iterations: 52

## a: 0.120285832190725  iterations: 52

## a: 0.12028358807766  iterations: 52

## a: 0.120349691629019  iterations: 53

а держится около 0.12 и количество итераций держится в одном интервале Количество итераций меняется при удалении значимых точек : чем ближе стартовое значение а к искомому, тем меньше количество итераций